Load all required libraries.

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages ---------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.0
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 3.6.3
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'readr' was built under R version 3.6.3
## Warning: package 'dplyr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## Warning: package 'plotly' was built under R version 3.6.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(broom)
## Warning: package 'broom' was built under R version 3.6.3

Read in raw data from RDS.

raw_data <- readRDS("./n1_n2_cleaned_cases.rds")

Make a few small modifications to names and data for visualizations.

final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
  rename(Facility = wrf) %>%
  mutate(Facility = recode(Facility, 
                           "NO" = "WRF A",
                           "MI" = "WRF B",
                           "CC" = "WRF C"))

Seperate the data by gene target to ease layering in the final plot

#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>% 
  select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke)) %>%
  group_by(date) %>% summarise_if(is.numeric, mean)

#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')
#remove facilty C for now
#only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
#only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]
#average N1 and N2
average_both <- only_positives %>%
  group_by(date, Facility) %>% 
  summarise(mean_copy_num_L = mean(mean_copy_num_L))
## `summarise()` regrouping output by 'date' (override with `.groups` argument)
        #renders the main plot layer three as positive target hits
        
        p77 <- plotly::plot_ly() %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = average_both,
                                       symbol = ~Facility,
                                       marker = list(color = '#E7298A', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
            layout(yaxis = list(title = "SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 type = "log",
                                 dtick = 1,
                                 automargin = TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #adds the limit of detection dashed line
        p77 <- p77 %>% plotly::add_segments(x = as.Date("2020-03-14"), 
                                          xend = ~max(date + 10), 
                                          y = 3571.429, yend = 3571.429,
                                          opacity = 0.35,
                                          line = list(color = "black", dash = "dash")) %>%
          layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y", 
                                    text = "Limit of Detection", showarrow = FALSE))

        

        p77
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Build the main plot

      #first layer is the background epidemic curve
        p1 <- only_background %>%
              plotly::plot_ly() %>%
              plotly::add_trace(x = ~date, y = ~new_cases_clarke, 
                                type = "bar", 
                                hoverinfo = "text",
                                text = ~paste('</br> Date: ', date,
                                                     '</br> Daily Cases: ', new_cases_clarke),
                                alpha = 0.5,
                                name = "Daily Reported Cases",
                                color = background_color,
                                colors = background_color,
                                showlegend = FALSE) %>%
            layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #renders the main plot layer two as seven day moving average
        p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke, 
                             type = "scatter",
                             mode = "lines",
                             hoverinfo = "text",
                            text = ~paste('</br> Date: ', date,
                                                     '</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
                             name = "Seven Day Moving Average Athens",
                             line = list(color = seven_day_ave_color),
                             showlegend = FALSE)
      

        
        #renders the main plot layer three as positive target hits
        
        p2 <- plotly::plot_ly() %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n1,
                                       symbol = ~Facility,
                                       marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n2,
                                       symbol = ~Facility,
                                       marker = list(color = '#D95F02', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
            layout(yaxis = list(title = "SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 type = "log",
                                 dtick = 1,
                                 automargin = TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #adds the limit of detection dashed line
        p2 <- p2 %>% plotly::add_segments(x = as.Date("2020-03-14"), 
                                          xend = ~max(date + 10), 
                                          y = 3571.429, yend = 3571.429,
                                          opacity = 0.35,
                                          line = list(color = "black", dash = "dash")) %>%
          layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y", 
                                    text = "Limit of Detection", showarrow = FALSE))

        

        p1
## Warning: Ignoring 1 observations
        p2

Combine the two main plot pieces as a subplot

p_combined <-
    plotly::subplot(p2,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
p_combined
p_combined_ave <-
    plotly::subplot(p77,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
p_combined_ave
save(p_combined_ave, file = "./plotly_fig_ave.rda")

Save the plot to pull into the index

save(p_combined, file = "./plotly_fig.rda")

Save an htmlwidget for website embedding

htmlwidgets::saveWidget(p_combined, "plotly_fig.html")

Build loess smoothing figures figures

#create smoothing data frames 
#n1
smooth_n1 <- only_n1 %>% select(-c(Facility)) %>% 
  group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
  summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
  ungroup() %>%
  mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
  mutate(target = "N1")
## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)
#n2
smooth_n2 <- only_n2 %>% select(-c(Facility)) %>% 
  group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
  summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
  ungroup() %>%
  mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
  mutate(target = "N2")
## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)
#add trendlines 
#extract data from geom_smooth
#n1 extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_n1 <- ggplot(smooth_n1, aes(x = date, y = log_sum_copies_L)) + 
  stat_smooth(aes(outfit=fit_n1<<-..y..), method = "loess", color = '#1B9E77', 
              span = 0.6, n = 142)
## Warning: Ignoring unknown aesthetics: outfit
#n2 extract
extract_n2 <- ggplot(smooth_n2, aes(x = date, y = log_sum_copies_L)) + 
  stat_smooth(aes(outfit=fit_n2<<-..y..), method = "loess", color = '#1B9E77', 
              span = 0.6, n = 142)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#n1
extract_n1
## `geom_smooth()` using formula 'y ~ x'

fit_n1
##   [1] 11.25620 11.38418 11.50944 11.63126 11.74894 11.86175 11.96899 12.06994
##   [9] 12.16507 12.25555 12.34168 12.42377 12.50215 12.57711 12.64898 12.71634
##  [17] 12.77787 12.83399 12.88510 12.93162 12.97397 13.01257 13.04784 13.08018
##  [25] 13.11001 13.13775 13.16383 13.18864 13.21261 13.22798 13.22961 13.22188
##  [33] 13.20914 13.19578 13.18616 13.18464 13.17997 13.15969 13.12641 13.08273
##  [41] 13.03126 12.97458 12.91530 12.85603 12.79936 12.74789 12.70423 12.67098
##  [49] 12.65072 12.64608 12.65183 12.66083 12.67280 12.68744 12.70446 12.72358
##  [57] 12.74450 12.77370 12.81645 12.87072 12.93445 13.00558 13.08207 13.16187
##  [65] 13.24292 13.32316 13.40055 13.47304 13.53858 13.59510 13.64056 13.68320
##  [73] 13.73164 13.78444 13.84016 13.89736 13.95458 14.01038 14.06332 14.11196
##  [81] 14.15484 14.19053 14.21758 14.23454 14.23997 14.23066 14.20729 14.17427
##  [89] 14.13601 14.09692 14.06140 14.03385 14.00450 13.96207 13.90865 13.84634
##  [97] 13.77727 13.70353 13.62724 13.55050 13.47541 13.40409 13.33865 13.28118
## [105] 13.23381 13.19863 13.16875 13.13700 13.10506 13.07463 13.04740 13.02505
## [113] 13.00927 12.99783 12.98754 12.97864 12.97136 12.96594 12.96263 12.96166
## [121] 12.96328 12.96772 12.97522 12.98603 13.00038 13.01851 13.04067 13.06594
## [129] 13.09343 13.12339 13.15608 13.19173 13.23061 13.27295 13.31876 13.36782
## [137] 13.42004 13.47532 13.53358 13.59472 13.65866 13.72530
#n2
extract_n2
## `geom_smooth()` using formula 'y ~ x'

fit_n2
##   [1] 11.01877 11.18829 11.35391 11.51501 11.67095 11.82109 11.96479 12.10143
##   [9] 12.23140 12.35574 12.47472 12.58862 12.69771 12.80227 12.90258 12.99738
##  [17] 13.08549 13.16727 13.24311 13.31338 13.37846 13.43871 13.49452 13.54625
##  [25] 13.59429 13.63901 13.68077 13.71997 13.75696 13.78560 13.80172 13.80874
##  [33] 13.81007 13.80911 13.80927 13.81395 13.81159 13.79030 13.75301 13.70265
##  [41] 13.64215 13.57444 13.50245 13.42910 13.35733 13.29007 13.23024 13.18078
##  [49] 13.14461 13.12467 13.11094 13.09317 13.07390 13.05567 13.04104 13.03255
##  [57] 13.03274 13.04510 13.07017 13.10605 13.15087 13.20274 13.25977 13.32008
##  [65] 13.38179 13.44301 13.50186 13.55645 13.60490 13.64532 13.67583 13.70316
##  [73] 13.73464 13.76932 13.80620 13.84434 13.88274 13.92044 13.95648 13.98987
##  [81] 14.01965 14.04484 14.06448 14.07759 14.08320 14.07835 14.06301 14.04059
##  [89] 14.01448 13.98809 13.96482 13.94809 13.93244 13.91074 13.88387 13.85270
##  [97] 13.81811 13.78098 13.74218 13.70260 13.66310 13.62456 13.58787 13.55389
## [105] 13.52351 13.49760 13.48084 13.47431 13.47349 13.47383 13.47080 13.45988
## [113] 13.43652 13.40569 13.37520 13.34474 13.31400 13.28269 13.25050 13.21714
## [121] 13.18229 13.14566 13.10695 13.06585 13.02206 12.97528 12.92521 12.87305
## [129] 12.81999 12.76568 12.70978 12.65195 12.59186 12.52916 12.46384 12.39621
## [137] 12.32634 12.25433 12.18028 12.10426 12.02637 11.94671
#assign fits to a vector
n1_trend <- fit_n1
n2_trend <- fit_n2

#extract y min and max for each
limits_n1 <- ggplot_build(extract_n1)$data
## `geom_smooth()` using formula 'y ~ x'
limits_n1 <- as.data.frame(limits_n1)
n1_ymin <- limits_n1$ymin
n1_ymax <- limits_n1$ymax

limits_n2 <- ggplot_build(extract_n2)$data
## `geom_smooth()` using formula 'y ~ x'
limits_n2 <- as.data.frame(limits_n2)
n2_ymin <- limits_n2$ymin
n2_ymax <- limits_n2$ymax

#reassign dataframes (just to be safe)
work_n1 <- smooth_n1
work_n2 <- smooth_n2

#fill in missing dates to smooth fits
work_n1 <- work_n1 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n1 <- work_n1$date
work_n2 <- work_n2 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n2 <- work_n2$date

#create a new smooth dataframe to layer
smooth_frame_n1 <- data.frame(date_vec_n1, n1_trend, n1_ymin, n1_ymax)
smooth_frame_n2 <- data.frame(date_vec_n2, n2_trend, n2_ymin, n2_ymax)
#make plotlys

#plot smooth frames
p3 <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_n1, y = ~n1_trend,
                    data = smooth_frame_n1,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n1,
                                  '</br> Median Log Copies: ', round(n1_trend, digits = 2),
                                  '</br> Target: N1'),
                    line = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
plotly::add_lines(x = ~date_vec_n2, y = ~n2_trend,
                  data = smooth_frame_n2,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n2,
                                  '</br> Median Log Copies: ', round(n2_trend, digits = 2),
                                  '</br> Target: N2'),
                    line = list(color = '#D95F02', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
plotly::add_ribbons(x ~date_vec_n1, ymin = ~n1_ymin, ymax = ~n1_ymax,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n1, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(n1_ymax, digits = 2),
                                  '</br> Min Log Copies: ', round(n1_ymin, digits = 2),
                                  '</br> Target: N1'),
                    name = "",
                    line = list(color = '#1B9E77')) %>%
plotly::add_ribbons(x ~date_vec_n2, ymin = ~n2_ymin, ymax = ~n2_ymax,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n2, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(n2_ymax, digits = 2),
                                  '</br> Min Log Copies: ', round(n2_ymin, digits = 2),
                                  '</br> Target: N2'),
                    name = "",
                    line = list(color = '#D95F02')) %>%
                layout(yaxis = list(title = "Total Log SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
    plotly::add_segments(x = as.Date("2020-06-24"), 
                                          xend = as.Date("2020-06-24"), 
                                          y = ~min(n1_ymin), yend = ~max(n1_ymax),
                                          opacity = 0.35,
                                          name = "Bars Repoen",
                                          hoverinfo = "text",
                                          text = "</br> Bars Reopen",
                                                 "</br> 2020-06-24",
                                          showlegend = FALSE,
                                          line = list(color = "black", dash = "dash")) %>%
    plotly::add_segments(x = as.Date("2020-07-09"), 
                                          xend = as.Date("2020-07-09"), 
                                          y = ~min(n1_ymin), yend = ~max(n1_ymax),
                                          opacity = 0.35,
                                          name = "Mask Mandate",
                                          hoverinfo = "text",
                                          text = "</br> Mask Mandate",
                                                 "</br> 2020-07-09",
                                          showlegend = FALSE,
                                          line = list(color = "black", dash = "dash")) %>%
    plotly::add_segments(x = as.Date("2020-08-20"), 
                                          xend = as.Date("2020-08-20"), 
                                          y = ~min(n1_ymin), yend = ~max(n1_ymax),
                                          opacity = 0.35,
                                          name = "</br> Classes Begin",
                                                 "</br> 2020-08-20",
                                          hoverinfo = "text",
                                          text = "Classes Begin",
                                          showlegend = FALSE,
                                          line = list(color = "black", dash = "dash")) %>%
      plotly::add_segments(x = as.Date("2020-10-03"), 
                                          xend = as.Date("2020-10-03"), 
                                          y = ~min(n1_ymin), yend = ~max(n1_ymax),
                                          opacity = 0.35,
                                          name = "</br> First Home Football Game",
                                                 "</br> 2020-10-03",
                                          hoverinfo = "text",
                                          text = "First Home Football Game",
                                          showlegend = FALSE,
                                          line = list(color = "black", dash = "dash")) %>%
  plotly::add_markers(x = ~date, y = ~log_sum_copies_L,
                      data = smooth_n1,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_sum_copies_L, digits = 2)),
                       marker = list(color = '#1B9E77', size = 6, opacity = 0.65)) %>%
    plotly::add_markers(x = ~date, y = ~log_sum_copies_L,
                      data = smooth_n2,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_sum_copies_L, digits = 2)),
                       marker = list(color = '#D95F02', size = 6, opacity = 0.65))

p3

Create final trend plot by stacking with epidemic curve

smooth_extract <-
    plotly::subplot(p3,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
smooth_extract
save(smooth_extract, file = "./smooth_extract.rda")